#include "field.h"

void field_laplacian(Field        *f_laplacian,
                     const Field  *f_source,
                     const double  inv_dx2,
                     const double  inv_dy2
                     )
{
    
    CHECK_FIELD(f_laplacian);
    CHECK_FIELD(f_source);
    CHECK_FIELDS(f_laplacian,f_source);
    {
        const int       nx        = f_laplacian->nx;
        const int       ny        = f_laplacian->ny;
        const int       nx1       = nx-1;
        const int       ny1       = ny-1;
        double **       laplacian = ARRAY_OF(f_laplacian);
        const double ** a         = CONST_ARRAY_OF(f_source);
        int i,j;
        
        for(j=0;j<ny;++j)
        {
            const int jm = (j <= 0)   ? ny1 : j-1;
            const int jp = (j >= ny1) ? 0   : j+1;
            
            for(i=0;i<nx;++i)
            {
                const int im = (i <= 0)   ? nx1 : i-1;
                const int ip = (i >= nx1) ? 0   : i+1;
                
                const double a0 = a[j][i];
                const double mid = -(a0+a0);
                
                laplacian[j][i] =
				inv_dx2 * ( a[j][im] + mid + a[j][ip] ) +
                inv_dy2 * ( a[jm][i] + mid + a[jp][i] );
            }
        }
    }
}
